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ABSTRACT 

It has been shown that the large-scale correlation functions of the density field (and velocity 
divergence field) follow a specific hierarchy in the quasilinear regime and for Gaussian initial 
conditions (Bernardeau 1992). The exact relationships between the cumulants of the probability 
distribution functions (the so-called S p parameters) are however sensitive to the smoothing window 
function applied to the fields. In this paper, I present a method to derive the whole series of the 
S p parameters when the density field is smoothed with a top-hat window function. The results are 
valid for any power spectrum and any cosmological parameters. Similar calculations are presented 
for the velocity divergence field. 

The resulting shapes of the one-point probability distribution functions of the cosmic density 
and the velocity divergence fields are given as a function of the power spectrum and £1. Simple 
analytical fits are proposed when the index of the power spectrum is —1. Comparisons with numerical 
simulations prove these analytical results to be extremely accurate. 



Subjects headings : Cosmology: theory - large scale structure of the universe - 
Galaxies - clustering 
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1. Introduction 

In most cosmological models, the large-scale structures of the Universe are thought to form 
by gravitational instabilities from small density fluctuations. In the mildly nonlinear regime, where 
the density fluctuations are still small, the mass distribution and motion of the Universe can reliably 
be described by the behavior of a pressureless fluid in an expanding Universe. It is then reasonable 
to expect perturbation theories to give a good description of the statistical properties of the large- 
scale cosmic fields. The linear theory has been widely evoked for the evolution of the rms density 
fluctuation and it has been also used for the velocity density relationship at large scales. However 
other large-scale features escape the first order perturbation theory, such as, for instance, the 
skewness of the density distribution function for Gaussian initial conditions. Indeed it is given 
by the second order perturbation theory (Peebles 1980). Actually, the characterization of the large- 
scale matter distribution and motion involve many more quantities than the simple rms density 
or velocity fluctuation. In the mildly nonlinear regime the shape the one-point density probability 
distribution function significantly departs from a Gaussian distribution. This evolution can be seen 
through the emergence of non-zero cumulants (of order greater than two) for this distribution, the 
skewness being the first one. 

The assumption of Gaussian initial conditions is crucial for the derivation of such quantities 
and in the following, I will assume that the initial fluctuations follow Gaussian statistics. The original 
derivation of the skewness by Peebles has been improved since to take into account the effects of 
filtering (that in practice cannot be avoided) with numerical (Goroff et al. 1986, Grinstein & Wise 
1987) or analytical (Juszkiewicz, Bouchet & Colombi 1993a, Bernardeau 1994b) calculations and 
the effects of redshift space distorsion (Bouchet et al. 1992) as well as its £1 and A dependences. On 
the other hand, Fry (1984), Fry, Melott & Shandarin (1993) explored the geometrical dependences 
of the three-point correlation function given by these perturbative calculations. Similar calculations 
for the divergence of the velocity field have also been done (Bernardeau et al. 1994, Bernardeau 
1994b) leading to new tools to constrain the cosmological parameters from the statistics of the 
large-scale cosmic velocity field. The kurtosis, obtained from the fourth moment of the distribution 
function, has been as well examined in details with perturbation theory. This quantity was originally 
derived by Fry (1984) and by Goroff et al. (1986) from third order perturbation theory. The effects 
of smoothing with a top hat window function on this quantity have been calculated analytically 
by Bernardeau (1994b) for both the density and the velocity divergence. For all these results the 
perturbation theory proved its efficiency when it has been compared with numerical simulations 
(previous reference and Juszkiewicz et al. 1993b). 

However, up to now, only the skewness and the kurtosis were known when the smoothing 
effects are taken into account, whereas Bernardeau (1992) was able to derive the full hierarchy 
of the large-scale cumulants when the smoothing effects are neglected. The knowledge of such 
a hierarchy makes possible the derivation of the shape of the one-point probability distribution 
function of the density or velocity divergence in a regime of low a. It is then of crucial interest to be 
able to derive the smoothing effects for the whole hierarchy of the cumulants. Other methods using 
various kinds of approximations have been used to derive the shape of this distribution function. 
Kofman (1991) and Kofman et al. (1993) use the Zel'dovich approximation (Zel'dovich 1970) for the 
large-scale dynamics to derive the shape of this distribution function but without taking into account 
the smoothing effects. The use of the Zel'dovich approximation, however, leads to wrong limits for 
the cumulants of the distribution as shown by Bernardeau & Kofman (1994) which undermines this 
method for a practical use. 

In this paper, I report how it is possible to obtain the full hierarchy of the cumulants when the 
smoothing effects (for a top-hat window function only) are taken into account. Part 2 is devoted to 
a recall of the results that have been obtained when the smoothing effects are neglected. In Part 3, 
I present the general method to get the smoothing corrections for the density field. The resulting 
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shape of the one-point probability distribution function of the density is examined in detail. Part 4 
is devoted to similar calculations for the divergence of the velocity field. 



2. The correlation hierarchy when the smoothing effects are neglected 

Most of the results I am going to present in this part have been derived in detail in a previous 
paper (Bernardeau 1992). The notation, however, will be slightly updated. 

As mentioned in the introduction, I assume that the large-scale structures of the Universe form 
by gravitational instabilities from Gaussian initial fluctuations and that the content of the universe 
behaves like a pressureless fluid. No assumptions are made about the cosmological parameters £1 
and A. I will also assume that the multi streaming regime that may appear at small scales does not 
change the large-scale mass distribution and motion. The motion equations are then supposed to 
be given by the single stream approximation, that is that, at any location in space one can associate 
a unique peculiar velocity. 

In the following the matter density field is described by the overdensity 6(x) = p(x)/p — 1 
and the peculiar velocity field, u(x), is, according to a general result of perturbation theory (Peebles 
1980), assumed to be irrotational and then can be described by its divergence, 

0(x) = -j-V.u(x). (1) 

2.1. Definition of the cumulants 

The basic objects of my analysis are the cumulants of the distribution function of 6 and 6. 
The p th cumulant of the distribution function of 6, (fi p ) c is defined recursively from the p th moment 
following the rules, 

<* 2 > e = <* 2 >^ 2 , 

= <* B >-io<* 2 > e <* 3 > e . 

and, in general, to obtain the p th cumulant one should consider all the decompositions of a set of 
p points in its subsets (but the one being only the set itself); multiply, for each decomposition, the 
cumulants corresponding to each subset and subtract the results of all these products thus obtained 
from the p th moment. It is worth noticing that p/p and 6 have the same cumulants. A similar 
definition applies for the cumulants of the distribution function of 6. The second moment, a 2 , of the 
density distribution is related to the power spectrum of the fluctuations, P(k). In the following I 
assume that the rms fluctuation is small enough to be accurately given by the linear theory so that, 

where D(i) is the time dependence of the growing mode and to the present time if P(k) is normalized 
to the present day large-scale fluctuations. The rms fluctuation of the divergence field, ag, is related 
to a through the function, 
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(a is the expansion factor) with 

<7 9 =/(fi,A)<7. (4) 

The function / can be accurately approximated by /(SI, A) = £1° 6 according to Peebles (1980) and 
Lahav et al. (1991). 

Each cumulant then contains a new piece of information on the statistical properties of the 
fields related to the shape of the distribution function. A Gaussian distribution, for instance, is 
characterized by only one non-zero cumulant, its variance (see Bernardeau & Kofman 1994 for a 
more detailed discussion on the relationship between the shape of the density probability distribution 
function and the cumulants in the mildly nonlinear regime). 

2.2. Large-scale behavior of the cumulants 

The starting point of this analysis are the motion equations that describe the behavior 
of a pressureless fluid in the single stream approximation, the continuity, Euler and Poisson 
equations. The density field 6(x) can be expanded with respect to the initial fluctuation field, 
(5(x) = X2p=i ^ P H X ) where 6^ is of order p with respect to the initial fluctuation field. Using the 
motion equation it is possible to determine 6^ order by order. For an Einstein-de Sitter universe 
the solution can be expressed in terms of the Fourier transforms of the initial fluctuations, <5(k) , so 
that (Goroff et al. 1986) 

^ 1)(X) = / (^^ k ) 6Xp ( lkx ) 

(5) 

6(P)(X) = J (^•••(^TI^ 1 )---^)^^ 1 '---'^) exp[i(k 1 + ... + k p )x] 

where #( p )(ki, . . . , k p ) is an homogeneous function of the wave vectors ki, . . . , k p . Then, the leading 
term of the cumulant (6 P ) c when a is small is given by 

<*>> e = £ <n^> c (6) 

combinations, p(i) 8 = 1 

where the sum is made over all the possible combinations p(i), i = 1 . . .p, for which 

p 

P(i)>l, £>(i) = 2p-2. (7) 

8 = 1 

This particular result is due to the hypothesis of Gaussian initial conditions. It is then possible to 
show that the terms in (6) are all products of the rms density fluctuation at the power 2p — 2 by 
some combination of the vertices, 

vM = <«Cp) )Jv*p = ^L-J d n 1 ...dn p #^)(k 1; . . . , k p ) (8) 

that is the angular average of the function #( p )(ki, . . . ,k p ). This angular average leads to keep 
only the monopole term of . The coefficients i> p are in general slightly time dependent but it 
just reflects a small change of these coefficients with the cosmological parameters (as shown in the 
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following), and this dependence vanishes for an Einstein-de Sitter Universe. Then, each cumulant, 
(fi p ) c , decreases with the variance like (6 2 ) P so that the ratios (6 P ) J (6 2 ) P have a finite limit, 

J7Z^T^S p (a) when (6 2 )^0. (9) 

The coefficients S p are closely related to the vertices v p . Technically it corresponds to a tree 
summation (see Fry 1984, Bernardeau 1992). It is possible to get a closed analytical relationship 
between the vertices and the cumulants at the level of their generating functions. Let me define 

rta,y) = y j S p (a) { -^—y p , (5 2 = 1), (10) 

7^2 p - 



and 



Then we have 



Gs(a,T) = ^Vpia) 



(-rf 



P =i 



<f(a, y) = yGs(a, t) - \yT^-Qs{a, r), r = -y^-Q 6 {a, r), 

Z QT QT 



(11) 



(12) 



These rules for tree summations have been given by De dominicis and Martin (1964) (see also, Des 
Cloiseaux and Jannink 1987) and has been used in an astrophysical context by Schaeffer (1985), by 
Bernardeau & Schaeffer (1992) and Bernardeau (1992) for the cumulants of the density distribution 
function. 

The last step of the derivation of the S p parameters is the derivation of the Gs{cl,t) function 
from the motion equations. Taking advantages of the geometrical averages in (8) it is possible to 
get a simple 2 nd order differential equation for Qs(a,r) that corresponds to the spherical collapse 
dynamics, i.e. , Qs is the density contrast of an object of linear overdensity — r. For an Einstein-de 
Sitter Universe the time dependence of the generating functions vanishes and the function Gs{j) is 
given by, when r < 0, 



9 (6 - sin ( 



2 (1 - cos( 



2/3 



(13a) 



and when r > 0, 



9 (sinh0 - 0) 2 _ 3 

Gs ~2 (cosh 9 - l) 3 ~ ' T ~5 



-i 2/3 



-(sinh ( 



(136) 



For SI ^ 1 and A = 0, the analytical solutions are given in Bernardeau (1992, eq. [35]). 
In general (A ^ 0) there is no analytical solution for the spherical collapse dynamics. However, 
whatever the cosmological parameters, the function Gs{j) is well described by the function, 

( 1 + y) ^ - 1 ' ( 14 ) 

showing that the time dependence of Gs{a,T) (and consequently the one of S p (a)) is extremely 
weak. In the following I will widely use the approximation (14). One may wonder why such an 
approximation is good since it even fails to reproduce accurately the position of the singularity. 
Actually in the system (12) the values of interest for r are r > —0.42 for which the approximation 
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is extremely good. This limit in r comes from the implicit equation which has a solution only in 
this domain, so that the vicinity of the singular point in (14) is not relevant. 

Similar properties can also be derived for the divergence of the velocity field. Thus the large- 
scale limit of the cumulants of the velocity divergence probability distribution function follows the 
rules 

(0P) 

N /c T p (a) when (6> 2 ) 0, (15) 



<* 2 > 



p-i 



and the generating function of T p (a), ipg(a,y) = Ylp^=2 ^p( a )f P ( — 1) P 1 /Pn i s ; as well, related to a 
generating function of vertices, 



Pe(a,y) = yGe(a,T) - ^yT-^-g (a,T), r 
z ar 



-y—6e(a,T), 
dr 



with, in general, 



g e [a,f(n,A)r] 



^&(a,r) + /(n,A)^&(a,r) 
aa ar 



/[l + &(a,r)]. 



As the time dependence of Qs(a, r) is weak we end up with, 

Ge[a, /(fi, A)r] = -f(il, A)^G s (a, r)/ [1 + G s (a, r)]. 

ar 

Using the approximative form of /(SI, A), we obtain 



(16a) 



(166) 



G s (a,T) = r 1 + 



2r 

m^ 1 



(17) 



2.3. Tie one-point probability distribution function 

Up to this point what have been obtained is an exact mathematical result for the large-scale 
limits of the cumulants. To construct the probability distribution function, I will now assume that 
the ratios {& p ) c / {& 2 ) P actually equal their limit for the scales of interest. This approximation is 
expected to be valid as long as a is small compared to unity since the first corrective terms for the 
ratios would be of the order of a 2 in a perturbative calculation. 

The construction of the one-point probability distribution function of the density or of the 
divergence from the generating function of the cumulants has been discussed at length and for various 
regimes by Balian & Schaeffer (1989). The general result is that, 



p(8)d8 



+loo 



dy 



O • 2 eX P 



<p(y) 

n-2 + rr2 



d6, 



(18) 



where the integration is made in the complex plane. The limit of interest for this analysis is for 
c — *■ 0. In such a case and provided 6 < 6 C (Bernardeau 1992, Eq. 46), 



p(6)d6 



1 



l-r^'(r)/^(r) 



27T<7 2 



1/2 



exp 



T 

"2^2 



dS, g s (r) = 6 ; 



(19) 



where 6 C ps 0.66 is the value of the density contrast that cancels 1 — tQ' s '(t)/Q' s (t). The function 
f{y) can as well be reobtained from the shape of the density probability distribution function, 



exp 



db p(8) exp 



yb_ 



(20) 
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The use of the form (18) for p(8) allows indeed to revover <p(y). But the approximate form (19), 
with the use of a second saddle point approximation, allows as well to recover the form (12) of 
f{y), showing that the form (19) still contains all the information on the generating function of the 
cumulants. This property will be of interest for the calculations in the next section. 

A similar derivation can be done for the shape of the one-point distribution function of the 
velocity divergence. 



3. The smoothing effects for the density field 



3.1. The case of the skewness 

To illustrate the nature of the problem I will consider in this paragraph the case of the 
skewness of the density distribution function. It has been studied in detail in various papers cited 
in the introduction. 



3.1.1. Calculation in Lagragian space 

Let me first start with a Lagrangian description of the dynamics. In such a case the local 
density is obtained by the inverse of the Jacobian of the transformation, J(q) = |9x/9q|, between 
the positions of the particles in Lagrangian space, q, and their positions in real space, x. The 
derivation of the third moment requires the calculation of the density field at the second order in 
the initial fluctuation field. The first order of the density fluctuation is simply given by — J^iq) 
and at the second order by <$( 2 ) = (J 1 - 1 - 1 ) — J 1 - 2 - 1 . We then have (Bouchet et al. 1992), 



* (2) (q) = |[* (1) (q)] 2 + fE^S(«0 

a, (3 



>S(q). 



(21) 



where <^(q) if the gravitational potential at the Lagrangian position q and the subscript |Cf| g means 
the spatial derivatives with respect to the components a and [3 of q. It is more convenient to write 
this expression in Fourier space. Let me denote <5(k) the Fourier components of the initial density 
fluctuations (normalized to the present time), then 



^ (1) (q) 
^ (2) (q) 



d 3 k 

(2tt) 3 / 2 

d 3 ki 



<5(k) exp (iq.k) 



d 3 k 2 



(2tt) 3 / 2 (2tt) 3 / : 



«(ki)«(k 2 ) 



ki.k 2 



(22) 



exp [iq.(ki + k 2 



The third moment of the density distribution is then given by 3(((5( 1 )) 6^) (which is its 
dominant term in perturbation theory). The calculation of this ensemble average reads (once the 
symmetry factors have been taken into account) 



<* 3 > = 6 



d 3 ki d 3 k 2 
(2tt) 3 (2tt) 3 



P(ki)P(k 2 ) 



ki.k 2 

k\k 2 



(23) 



The integral over the angle between ki and k 2 in (23) is simple in the absence of window function 
and leads to the result, 



34 



<* 3 > = y<* 2 > 



(24) 



3.1.1. Filtering in Lagragian space 
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Let me then consider the skewness of the density distribution function for a filtering at a 
given mass scale. The density is simply defined as the ratio of the initial volume by the final volume 
occupied by a certain amount of matter. I will only consider the case of a top-hat filter that means 
that I will examine the statistics the inverse of the volume of particles that lie, in Lagrangian space, 
in a sphere of given radius Ro- The volume occupied by these particles all along the evolution is 
given by the integral of the Jacobian over the sphere of radius Ro- The top-hat filter should then 
be applied to the field J(q) and not to the field 6(q). Once it has been done the filtered density 
field at first and second order read 



4%) 



4 2 } (q) 



: 8(k) exp (iq.k) W TR (kR ) 



d 3 k 

(2tt) 3 / : 

d 3 ki d 3 k2 
(2tt) 3 / 2 (2tt) 3 / : 



■«(ki)«(k 2 )exp[iq.(ki +k 2 



x W^ih R )W TE (k 2 R ) 



7 



1 



ki.k 2 

k x k 2 



W XH (|ki +k 2 | R ) 



(25) 



where W X jj(fc Ro) is the Fourier transform of a top-hat window function, 

3 



W TE (kR ) 



[sin(fci?o) — kRg cos(kRo)] 



(kR ) 3 

The integral over the wave vectors ki and k 2 is dramatically changed in, 



v4„> = 6 



d 3 ki d 3 k 2 



(2tt) 3 (2tt) 



3 P(*i)P(* 2 )W XH (*i R )W TE (k 2 R ) 



W TR (h Ro)W TR (h Ro) ~ ^W TE (|ki + k 2 | Ro) 



k\k 2 



(26) 



(27) 



A geometrical property of the top-hat window function given by Bernardeau (1994b, Eq. [A. 5]), 
however, shows that the integral over the angle between the wave vectors is also simple and leads 
exactly to the same coefficient as previously, 



showing that the filtering at a given mass scale does not change the skewness. 



3.1.3. Filtering in Eulerian space 



In practice, however, we are more interested in the behavior of the distribution function in 
Eulerian space, for a given smoothing radius, Ro. When the density field is expressed in term of the 
Eulerian coordinates, x, rather than the Lagrangian coordinates, its first and second order read, 



£W(x) 
^ 2 )(x) 



d 3 k 

(2tt) 3 / 2 

d 3 ki 



<5(k) exp (ix.k) 



d 3 k 2 



(2tt) 3 / 2 (2tt) 3 / 2 



<5(ki)<5(k 2 )exp [ix.(ki + k 2 



ki-k 2 

k\ 



ki.k 2 

k\k 2 



(28) 



The new term in 6^ comes from the expression of ^-^(q) and is due to the change of variable from 
q to x. If one calculates the skewness from the expression (28), that is without smoothing, one gets 
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the coefficient 34/7 as in equation (27) for the Lagrangian case. The smoothing of the density field 
changes, however, the result. Indeed, after filtering, the third moment for the density field is given 

by, 

<^o> = 6 / ^^P(h)P(k2)W TR (k 1 R )W TR (k 2 R ) 



x W XH (|ki + k 2 | R ) 



5 k!.k 2 , 2 /k!.k 2 x 2 ' 



7 kf 7 V hk 



(29) 



This expression can then be calculated using a second property of the top-hat window function 
given by Bernardeau (1994b, Eq. [A. 6]). It leads to 

<^>=(y+7i)<^> 2 , (30«) 



with 



7. = (*»> 

d log R 



There is here a corrective term that depends on the scale dependence of rms fluctuation. 
The origin of this correction can now be understood as a shift in mass scale. Indeed, for a given 
mass scale the skewness remains uncorrected as the analysis in the Lagrangian space demonstrated 
it. However, the filtering process in Eulerian space mixes different mass scales: a given density p 
measured at a radius Ro corresponds to a mass scale of p 4itRq/3 (and not 47ri?o/3). When the 
rms fluctuation decreases with mass, less high density regions and more low density regions can be 
expected than if a were scale independent. The skewness being sensitive to the asymmetry between 
the large density contributions and the low density contributions is then expected to decrease. The 
result (30) quantifies this departure. 

3.2. The correlation hierarchy with the smoothing effects 

3.2.1 Filtering in Lagrangian space 

Once again I will start by the density distribution in Lagrangian space, that is for a given 
mass scale. In a Lagrangian description the density is given by the inverse of the Jacobian. One can 
write any order p of the expansion of the Jacobian with the initial displacement field, 

j(P>it) = J ^^■■■■^^■■■^J^U---^)^pM^ + --- + K)] 

where j( p )(ki, . . . , k p ) is a function of the angular parts of the wave vectors ki,...,k p only. 
For instance for p = 2 we have j( 2 ) = 2/7(1 — u\ 2 ) with ui 2 = ki.k 2 /(&i k 2 ). To get the 

(p) 

expression of J^ o when the field is smoothed at the radius Ro one has to multiply the function 
j( p )(ki, . . . , k p ) by the Fourier transform of the window function applied to the sum of the wave 
vectors, VF X jj [|ki + . . . + k p | Ro]. The result obtained in the appendix C of Bernardeau (1994a) is 
a geometrical property of the function j( p )(ki, . . . ,k p ). Indeed the integral over the angular parts 
of the wave vectors gives 

/ dfii...dfi p W XH [|ki + ... + k p | R ]j (p \k 1 ,...,k p ) = 

1 (t) (31) 
( 47r)P ^w XH (ftiiJo)...W XH (ft p iJ ). 
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where j p contains only the monopole term of j( p )(ki, . . . , k p ), i.e. , the quantity you would have 
obtained simply by dropping the function W^jj in (27). The integral over the angles then lead to a 
factorization of the k\, . . . , k p dependence and the result is exactly what would have been obtained if 
the filtering had been made on the initial field that is replacing P(k) by W 2 (k) P(k). The continuity 
equation, together with the properties of the generating functions given by Bernardeau (1992, Eq. 
[26]), gives the generating function of the vertices for the density field, 

^W=(l + E^ i 7r) - 1 (32) 

which is exactly the function (11) obtained when the smoothing was neglected. So we show here that 
the filtering preserves the shape of the high-order correlations and does not change the generating 
function of the vertices. The property of the previous paragraph can then be generalized: the 
smoothing does not change the values of the S p parameters when it is done at a given mass scale. 



3.2.2. Filtering in Eulerian space 

The difficulty then lies in the derivation of the S p parameters for the smoothed density field 
in real space. From the analysis of the skewness it appears that it will be due to a mixing of 
different masses. In general the mapping between the Lagragian space and the Eulerian space is 
rather complicated. However, an exact relationship between the probability distribution functions 
of the filtered density field in Lagrangian and in Eulerian space is not necessary. What is needed is 
a trick to get the generating function of the large-scale cumulants that takes into account the mass 
mixing. Let me choose at random a point of matter xo in the field. I can define Pec(Ro, M)dM the 
probability distribution function of the mass contained within the radius Ro centered on xo. On the 
other hand, for a given mass Mo it exists a unique radius R so that the mass contained in a cell of 
radius R centered on xo is exactly Mo- Let Plc(Mo, R)dR be the probability distribution function 
of such a radius. Then the probability that a given volume Vo of radius Ro, contains an amount of 
matter greater than Mo is also the probability that the mass Mo occupies a volume smaller than 
Vo- As a result we have 

rco /■Ro 

/ p Ec (R ,M)dM= p Lc (M ,R)dR. (33a) 

Both pec and p^ c can be seen as density probability distribution functions for respectively a fixed 
radius Ro or a fixed mass Mo, so that the previous relation also reads 

/>oo />oo 

/ PEc (Ro,6)d6= / PLc (Mo,6)dS, (336) 

J 6 J 6o 

with So = 3Mo/(47ri?op) — 1- For a continuous field the relationship between the function pec and the 
Eulerian smoothed density distribution function ps is simple and reads, pec(^)^S = (1 + S) Pe(^)^S. 
In general the relation between p^ c and the Lagrangian smoothed density distribution function p^ 
is complicated. This is due to the fact that the particles that were initially in a spherical region end 
not necessarily in a spherical region in Eulerian space. However, as recalled in §2 the S p parameters 
are only determined by the spherical collapse dynamics, and as seen for the skewness the smoothing 
corrections come from the mass scale mixing. The possible aspherical movement of the matter is 
thus not expected to change the value of the large scale limit of the S p parameters. I will then 
neglect it and assume that p^ c = p^. Using the relation (33b) it is then possible to relate p E and 
PL- By differentiating with respect to the threshold So for a fixed radius Ro, we get 

f°° d 

(1 + 6o)pe(6o) = Pl(6o) + / pVo— PL (S)dS. 

JSo dM 
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The large-scale values of the S p parameters of pL(8)d8 are known and are the ones of the unsmoothed 
density field. We can then use the relation (18) to derive the expression of the density distribution 
function in Lagrangian space. We then obtain for the distribution function in Eulerian space, 



Pe(6o) 



1 



x exp 



yS 



da 2 (Mq) 
dM 



(34) 



Let me replace the mass dependence of a by a scale dependence, so that 3p Voder 2 (Mo )/dMo = 
d<7 2 (i?o)/dlog(i?o)- The generating function of the cumulants for the smoothed density field at a 
given scale Rq can then be obtained directly from the equation (34), 



exp 



dSps^S) exp 



by 



<r 2 (i?o)J 



(35) 



Actually we do not need to calculate exactly the expression of the right side of (35) but only the limit 
of its logarithm when a = 0. One can then use safely the saddle point approximation to integrate 
over y to get pe(J>) from (34). It leads to 



exp 



Ro da 2 (R ) 



<7 2 (i? )J 3(1 + 6 ) dR 
l-rQ' s '(r)/Q' s (r)^' 2 



_2Tra 2 ([l + 6] 1 / 3 R ) 
with Gs(t) = 6. 



exp 



6 + 
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Y 



x -d6 

i WP) 

6y t 2 

L"<7 2 (i? ) 2,7 2 ([l + ( 5]l/3 jRo ) 



Let me then define the function Qf by the implicit equation, 



G1(t s ) = Qs 



a(\l + Q S s (r s )]^R ) 
<t(Ro) 



(36) 
(37) 



(38) 



When t s = t a(Ro)/c ([1 + Qf (r <s )] 1 / 3 i?o) the term under the exponential reads [6y + 
(t s ) 2 / c 2\ / a 2 (Rq) , Qf (t s ') = 6. The integration over 6 of (37) can be made by a second saddle 
point approximation (similar to the one used to obtain the expression (11) from the approximate 
form (19)). The saddle point is then given by, 



(39) 



where t s is solution of the implicit equation, 



■ dr s 

~db 



(40) 



Once the logarithm has been taken, the limit a = conserves only the term that is under the 
exponential at the saddle point position so that <Pe(u) = yGf (t 6 ') + ( r<S ) 2 /2 and t s is given by 
the saddle point position (40). As a result the series of the parameters 5^, when the smoothing 
corrections are taken into account, is still obtained by a generating function of vertices. The function 
Qs has simply been replaced by Qf defined by the equation (38). 



power law index n 



power law index n 



Fig 1. The coefficients Sp (left panel) and \Tp\ (right panel) (Eqs. [9,15]) 



as a function of the power law index 



for p = 3, . . . , 7 (from bottom to top). The dashed lines correspond to the coefficients 5*3 and T3 obtained for a 
Gaussian window function by respectively Juszkiewicz et al. (1993a) and Bernardeau et al. (1994). 



The generating function of the cumulants of the smoothed density field is thus given by the 
relation, 



In the following I use the superscript s for the quantities that have been obtained when the filtering 
effects are taken into account. The relation (41) is the central result of this paper. All the applications 
that are presented in the following derive from the equations (41), (38) and (13). 

3.3. Example of cumulants 

By expanding the function f S {y) with respect to y it is then possible to recover the full 
hierarchy of the cumulants. We have already noticed that they will be changed due to the variation 
of a with scale. If a were scale independent, which corresponds of a power law index n = —3, then 
the cumulants are identical to the ones found without smoothing. In this paragraph I will express 
the dependence of a with scale through its logarithmic derivatives with scale, 



ip s {y) = yQ s d {T)-\yT^-g s d (T), r 
z ar 



(41) 



1 P = 



dlog p R 



(42) 



A tedious calculation (to be done with a symbolic calculator), can then give the expression 
of the first few S p parameters as a function of j p , 
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S3 = y +7i, 

= 60712 62 7l 7 7l 2 2 72 

4 1323 3 3 3' 

_ 200575880 1847200 71 6940 7 i 2 235 7 i 3 1490 72 50 7 i 72 10 73 

5 ~ 305613 + 3969 + 63 + 27 + 63 + 9 + 27 ' 
5 6 = 12650 + 12330 7l + 4512 7i 2 + 734.0 71 3 + 44.81 71 4 + 775.8 72 

+ 375.9 7i 72 + 45.56 71 2 72 + 3.889 72 2 + 20.05 73 + 4.815 71 73 + 0.1852 74, 

SV = 307810 + 383000 71 + 190700 71 2 + 47460 71 3 + 5914 71 4 + 294.8 71 5 

+ 27340 72 + 20300 71 72 + 5026 71 2 72 + 414.8 71 3 72 + 358.1 72 2 

+ 88.15 71 72 2 + 902.6 73 + 443.3 71 73 + 54.44 71 2 73 + 7.778 72 73 + 14.20 74 

+ 3.4577174 + O.O864275, 

These results are valid for an Einstein-de Sitter Universe and have been obtained from the expression 
(13) from ^(r). 

The first two coefficients have already been given by Bernardeau (1994b) from a direct 
calculation using perturbation theory up to the third order. It confirms that the method presented 
here, despite its technical simplicity (and the lack of a rigorous derivation), can give straightforwardly 
the full series of the coefficients. These coefficients are presented as a function of n assuming a power 
spectrum of index n in Fig. 1. Note that when n is close to 1 they vanish. In such a regime the 
results presented here should be used with caution: the high-order cumulants are expected to be 
dominated by higher order contributions of the perturbation theory. 

3.4. The resulting shape of the density distribution function 




Fig 2. The one— point probability distribution function of the density contrast in the quasilinear regime. In the left 
panel (7 = 0.5 and the power law index varies, n =0 (dotted line), -1 (dashed line), -2 (long dashed line), -3 (solid 
line). In the right panel, n = — 1 and (7 varies (<7 =0.3, 0.5, 0.8). In both cases the thin dashed lines correspond to 
the lognormal distribution (Eq. [52]). 
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Let me first make a remark on the theoretical expression of the density distribution function 
is Eulerian space. For the calculation of pe(J>) from Pl{6) in (34) it was assumed that the ratios 
(^ P ) c /(^ 2 ) 2 are cons tant in Lagrangian space. This form, however, does not imply that the similar 
ratios obtained for the density distribution function in Eulerian space follow the same scaling rules, 
although their limits when a = are the exact ones. In the absence of any other physical arguments 
I will assume in the following that these ratios are actually constant in Eulerian space to derive 
the density probability distribution function. That means that the expression (18) (with tp s (y) 
instead of f{y)), rather than the expression (34), should be used to derive the shape of the density 
probability distribution function. Comparisons with numerical simulations presented in the next 
subsection prove that the hypothesis of constant ratios is extremely powerful. 

For the sake of simplicity, I assume in the following that the scale dependence of the rms 
fluctuation can be approximated by a power law behavior, a(Ro) oc R ( n + 3 )/ 2 j n sucn a case the 
power spectrum would have a power law behavior, 

P(k)xk n . (43) 

This is a useful approximation which is accurate for most of the power spectrum of interest. (It is 
however possible to do the subsequent calculations without it.) In such a case, ji = — (n + 3), j p = 
when p > 2, and so 

gf(r) = g s (r[l + gf(r)]- (n+3)/6 ). (44) 
The use of the approximation (14) for the function Gs{j) leads to 

r=l(l + gff +3,,6 [(l + gfy 2,3 -l}. (45) 

Note that for n = —1 this relation can be inverted in Qf = [— r/3 + (1 + r 2 /9) 1 / 2 ] 3 — 1. Integrating 
the equation (18) one can obtain the shape of the probability distribution function for different 
values of n and of a. Some of these functions are given in Fig. 2. 

It is also possible to find analytical approximations of the density probability distribution 
function. One key feature to get these forms is to determine whether the equation rQf {r)/Qf (r) = 
1 has a solution or not. As can be checked from equation (45) this equation has a solution only 
if n < 0. Note that this limit actually is not due to the approximation (14) but is exact for all 
cosmological parameters. In table 1, I present the value of r that is solution of this equation and 
the values of the parameters of interest at this point. It corresponds indeed to a singular point in 
y = y s for the function tp s (y) and close to this value <p(y) behaves like, 

V s (y)- Vs ~-a s (y-y s fl\ (46) 

The values of y s , a s and <p s are given in table 1 for different values of n and will be of interest in 
the following. 

Table 1. Tie parameters of the critical point as a function of the spectral index, n, for the density 
distribution. 



n 


8c 


ys 


(Is 


<Ps 


-3 


.656 


-0.184 


1.84 


-0.030 


-2.5 


0.804 


-0.213 


2.21 


-0.041 


-2 


1.034 


-0.253 


2.81 


-0.058 


-1.5 


1.443 


-0.310 


3.93 


-0.093 


-1 


2.344 


-0.401 


6.68 


-0.172 


-0.5 


5.632 


-0.574 


18.94 


-0.434 
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Fig 3. The one— point probability distribution function of the density contrast. The thick solid line has been 
obtained from the equation (18), the long dashed lines correspond to the approximations (47) and (50) and the short 
dashed line to the approximation (48). From top to bottom and left to right the values of n are successively -3, -2, -1 
and and (7 = 0.5 in all cases. 

There are mainly two different approximations for the density distribution function depending 
on the value of 6. When 6 is lower than the critical value b c (that depends on n), one can write an 
expression similar to equation (19), 

exp(-^)d«, gf(r) = 6. (47) 



P f(8)d8 



1-rgf (r)/gf'(r) 

27T<7 2 



When the density is greater the 6 C another approximation can be used that is valid in the large 
density limit. The behavior of the density distribution in such a case in driven by the shape of the 
function f S {y) around its singularity and it reads (Balian & Schaeffer 1989, Bernardeau 1992, Eq. 
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[49]), 



P s 2 (8)d8 



3a s cr 



(l + ,5)- 5 / 2 exp [-\y s \8/a 2 + \p s \/a 2 ] d6. 



(48) 



When n > the expression (47) is valid whatever the value of 8, otherwise the expression (48) shows 
that the large density cut-off is in exp(— \y s \8/a 2 ) which is very different from a Gaussian like cut-off. 

In Fig. 3, I present these approximations compared to an exact integration of the equation 
(18). Note that the algorithm for the integration in the complex plane roughly follows the same 
kind of distinction for large and small values of 8. The basic idea for such an algorithm is to define 
a path for y in the complex plane so that <p(y) — 8y is always a positive real number thus assuring a 
rapid cutoff of the integral in (18). The value of y to be chosen for the intersection of the path and 
the real axis is either given by the saddle point y 
singular point y s when 8 > S c . 



-r/Gf '(t), gf(r) = 8 when 8 < 8 C or by the 



General features can then be presented at this stage. From the form (47) one can infer the 
form of the cut-off at low density. We obtain, 



p s (8)d8 



1 



4(27T<7 2 ) 1 /2 



2(4- n) 



1 



1/2 



+ 



exp 



9 (1_ 
"8 



8) 



d8. 



(49) 



It shows that the cut-off is sharper when n is small. On the other hand, table 1 shows that the 
cut-off is more gentle in the high density limit when n is small. In other words, the distribution is 
more Gaussian like when n is large. This is not a surprising result knowing the dependence of the 
S p parameters that are all decreasing functions of n (see Fig. 1). 

As can be seen in Fig. 3, the form (47) is extremely accurate for 8 ^ 8 C whereas the form 
(48) is not very accurate for large densities. It is then of interest to estimate the first corrective 
term of this form when the large density limit is released. Actually such a term would contain an 
extra factor 8~ 1 I 2 compared to the dominant term. I then propose an approximation based on the 
form (48) that is more accurate for large values of 8. In such a case, the density distribution can be 
written, 

pf (8) = [1 + 2(0.8- a)a- 13 p- 05 ]p s 2 (8). (50) 



The corrective factor has been constructed empirically from the shape of the probability 
distribution function calculated numerically from the form (18). It is quite accurate as long as 
a ^ 1. It is possible to get a reasonably good analytical fit of the shape of the distribution function 
by an interpolation between the form (47) and the form (50). for instance we can use 



p s (8) tap* (6) + 1(6) [pf(«)exp [0.05(1 + 5)" 



1(8) 



exp(p) 



exp(p) + exp(—p) 



P=(6-8 ) 



CO 



(51) 



In the form (51), the factor exp[0.05(l + 8)~ 4 ] has been put just to regularize this form for small 
values of 8. The values of So and <7o can be adjusted to get a smooth interpolation between the two 
forms. For instance we can take So = 0.2, 0.4, 2, <7o = 2, 1.5, 0.8 for respectively n = —3, —2, —1. 

Fig. 2 also shows one striking feature for the n = — 1 case. It turns out that for values of 
interest for the density and when a is lower than unity, the probability distribution function is quite 
close to a lognormal distribution 



Pi(6)d8 



1 



(2^) 1/2 



exp 



(ln[l + 5] + ^ 1 2 /2) / 
2<r? 



d8 
1+7' 



ln(l + ,r 2 ). 



(52) 
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This distribution was proposed by Coles & Jones (1991) as a good empirical fit for the density 
probability distribution function. This form, however, does not exhibit the same cut-offs both at 
large densities and small densities as the distribution functions obtained with the form (18) (Eqs. 
[48, 49]). That the lognormal distribution is close to the exact behavior when n = — 1 seems to be 
a mere coincidence (see Bernardeau & Kofman 1994 for more detailed discussions on this subject). 

3.5. Comparisons with numerical simulations 



p(<5) 




-10 12 3 

6 



P(<5) 




Fig 4. Comparison with numerical simulations for the one— point probability distribution function of the density 
contrast. The upper panel corresponds to the timestep a/ao = 0.6 and the lower panel to the timestep a/ao = 1 of 
a numerical simulation with CDM initial conditions. The density probability distribution function has been measured 
at three different radius, R = 5 h ^Mpc (circles), R = 10 h ^Mpc (squares), R = 15 h ^Mpc (triangles) in each 
case. The properties of the local density fluctuations, n and (J, are given in table 2 for each theoretical curve. 

Juszkiewicz et al. (1993b) performed a series of numerical simulations to test predictions of 
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perturbation calculations both on the density field and on the velocity divergence. The filtering, 
however, was made by a Gaussian window function instead of a top-hat window function which is 
likely to change the results. In Fig. 1, I present the theoretical value of S3 obtained by Juszkiewicz 
et al. (1993a) in case of a Gaussian window function. When n = —3 the filtering does not introduce 
any correction either for a top-hat or for a Gaussian window function and the smoothing effects are 
still quite comparable when n is lower than -2. But for values of n or the order of —1 or greater, the 
smoothing corrections are less important for a Gaussian window function than for a top-hat. To 
compare the analytical results with numerical simulations one can then find a sort of effective index 
which would give the correct skewness with a top-hat filter. For instance, for n = — 1 the effective 
index would be n e ff —1.5, for n = 0, n e ff. ~ —1.2. This recipe allows to estimate the value of 
S4 for a Gaussian window function. For n = —1 and n = we obtain respectively S4 20.1 and 
S4 16.2 in good agreement with the result obtained by Catelan & Moscardini (1993) with Monte 
Carlo numerical integrations. The theoretical probability distribution functions calculated with the 
effective indices prove to be extremely accurate to predict the shape of the measured probability 
density distribution functions presented in the paper of Juszkiewicz et al. (1993b). 

I also took advantage of a simulation with CDM initial conditions, kindly provided by H. 
Couchman (Couchman 1991), to compare the theoretical predictions with a realistic model. The 
simulation contains 2.1 10 6 particles, has CDM initial conditions, Ho = 50 km s _1 Mpc -1 , SI = 1 and 
was performed in a cubic box of 200 Mpc side with an adaptive P 3 M code. The measurements 
of the shape of the probability distribution function of the density contrast have been made at two 
different timesteps corresponding to an expansion factor of a/ao ~ 0.6 and a/ao = 1 where do is the 
final expansion factor (at the end of the simulation we have ag = 0.97). 

The values of the local indices n can be determined from the shape of the initial spectrum for 
each given smoothing radius. I then approximate the scale dependence of the rms density fluctuation 
at each smoothing radius by a power law given by the local index. For each timestep I measured 
the counts in cells probabilities on a grid of 50 3 points for three different cell radius. For such radius 
greater then 5 /i _1 Mpc the shot noise is negligeable in this simulation so that the counts in cells give 
the shape of the numerical density probability distribution function. The rms density fluctuation 
a is measured in each case. The values obtained for n and a are given in table 2. Knowing n and 
it, it is then possible to calculate the expected shape of the probability distribution function of the 
density contrast and to compare it to the numerical results. These comparisons are presented in 
Fig. 4. They show a perfect agreement between the theoretical predictions and the measurements 
in the simulations. The error bars have been obtained by dividing the simulations in 8 parts and by 
calculating the rms fluctuation between the different subsets for each smoothing radius and density 
contrast. These errors are dominated by finite size effects and are thus certainly overestimated for 
the measurements in the total sample. Perturbation theory proves to be extremely accurate even 
when 8 is by far larger than unity, and for values of a slightly above unity. 

4. The smoothing effects for the velocity field 

In this part I consider the properties of the probability distribution function of the large-scale 
velocity divergence field, #(x) (Eq. [1]). 

4.1. The generating function of the cumulants with the smoothing effects 

Let me comment a bit more on the analytical derivations presented in the previous part for 
the density field. The central argument is that the function Qs obtained without taking into account 
the smoothing effects is exact for a given mass scale. Then the resulting expression of the function 
Qf (r) is simply given by the function Qs{r) but calculated at the mass, Mo, of interest. The relevant 
value of r depends on the mass scale through the function a in a way that can be deduced from 
the equation (15) by identifying the terms under the exponential in both Lagrangian and Eulerian 
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case (which gives t s /c(Mo) = r/cr(Ro)). The relationship between the mass Mo and the scale Ro 
is directly related to the value of Gf(r s ) (see again Eq. [15] that gives 1 + Gf (t) = 3M /(47ri? 3 )). 
Then we obtain the equation (38). 

A similar analysis can be done for the cumulants of the divergence of the velocity field. They 
are not affected by smoothing in Lagrangian space for a given mass scale (this result can be obtained 
in a similar way than for the density field since the expression of the velocity divergence as a function 
of the displacement field involves the same geometrical expressions). The generating function of the 
vertices for the divergence can then simply be calculated at the mass scale of interest, that is at the 
same point r as previously, 



gf(f(ci,A)r s ) = g e 



(53) 



<t(Ro) 

With (18) and (45) it gives the expression of Qf as a function of r. 

4.2. Examples of cumulants 

The expression of the generating function of the coefficients T p can naturally be derived from 
the expression of Qf , 

<Pe (a, 2/) = yQe ( a > r ) ~ \ yT ~^^ ( a ' r )' T = -Vj^Sf (a, t). (54) 

The expansion of this function with respect to y then can give the expressions of the coefficient 
T p as a function of the logarithmic derivatives of a. I give here their expression for an Einstein-de 
Sitter Universe, 

26 

-T 3 = — + 7i , 



T 4 



5 



12088 338 7 i 7 7 i 2 2 72 

441 + 21 + 3 + ~3~ ' 
94262120 161440 7 i 260 71 2 235 71 3 130 72 50 71 72 10 73 

305613 + 567 + 3 + 27 + 7 + 9 + 27 ' 



T 6 = 4694 + 5951 71 + 2802 71 2 + 581.2 71 3 + 44.81 71 4 + 480.3 72 

+ 297.1 71 72 + 45.56 71 2 72 + 3.889 72 2 + 15.82 73 + 4.815 71 73 + 0.1852 74, 
-T 7 = 90310 + 146000 71 + 93680 71 2 + 29800 71 3 + 4704 71 4 + 294.8 71 5 

+ 13380 72 + 12710 71 72 + 3991 71 2 72 + 414.8 71 3 72 + 284.1 72 2 + 88.15 71 7 2 2 

+ 562.9 73 + 351.5 71 73 + 54.44 71 2 73 + 7.778 72 73 + 11.23 74 + 3.457 71 74 + 0.08642 75, 



These functions are presented in Fig. 1 for power law spectra. Note that the first two 
have already been derived with direct perturbative calculations (Bernardeau 1994b). In general 
the dependence of the coefficients T p with the cosmological parameters is simple and can be easily 
obtained from the dependence of Qg with the cosmological parameters. We then approximately have 

T p (Q,A)k 1 T P (Q = 1,A=0). (55) 

The use of the SI dependence exhibited by these coefficients has already been proposed to determine 
SI from the large-scale velocity fields (Bernardeau et al. 1994, Bernardeau 1994b). This dependence 
will also give a strong SI dependence for the shape of the probability distribution function. 
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Fig 5. The one— point probability distribution function of the divergence of the velocity field for n = — 3 (left panel) 
and n = — 2 (right panel). The thick solid line is the exact result of the quasilinear regime for an Einstein— de Sitter 
Universe. The short dashed lines correspond to the approximations (51) and (58). 

4.3. The shape of the probability distribution function 

A general rule can be given for the dependence of the probability distribution function with 
the cosmological parameters. It takes advantage of the dependence of the cumulants with SI and A 
(Eq. [50]) that leads to, 

AO 



:(f(n,A),9,a s )d9 = p* 



1, 



(?6 



/(n,A)'/(n,A) 



/(n,A)' 



(56) 



Table 2. The parameters of the critical point as a function of the spectral index, n, for the distribution of the 
velocity divergence. 



n 


Oc 


Vs9 


a s 


•t>s9 


-3 


-0.750 


-0.222 


1.67 


-0.042 


-2.5 


-0.961 


-0.266 


1.99 


-0.061 


-2 


-1.348 


-0.332 


2.51 


-0.100 


-1.5 


-2.325 


-0.447 


3.63 


-0.201 



The probability distribution function can then be calculated exactly the same way as the one 
of the density contrast, Fig 5. Once again two different analytical approximations can be derived. 
One is valid when the divergence is greater than the critical value 6 C (depending on n, see table 2) 
and it is given by 



_J l-rgf»/gf'(r) 

-gf'(r) i™\ 



Pf e (6)d6 



1/2 



exp 



2<t2 



de, gf(r) = e. 



(57) 
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The other one is valid when the divergence is lower than the critical value, 

P2e(0)d0 = ^^(3/2 " 0)- 5 ' 2 exp [\y, e \6/a 2 e + M/* 2 e ] d9, (58) 

where the parameters a s g, y s g and <p s g are given in table 2 for few values of n. 

It appears that when n > — 1 there is no singular point for the function <pf (y) so that the 
expression (57) is always valid. 

The form (57) is extremely accurate for large values of 6 whereas the form (58) underestimates 
the exact derivation of the probability distribution function. I then propose a correction to the form 
(58) which is in a better agreement with the exact form, 

pi e (0)de = [i + 30(.s - <T S )a; 13 (i.5 - ey 15 ] P s 2e (e)de. (59) 

In Fig. 5, the shape of these approximative forms are shown. Note that these forms have been 
given for an Einstein-de Sitter Universe but can be easily extended to any cosmological model with 
equation (56). 

The case n = — 1 is worth investigating for two reasons. First it corresponds roughly to 
the observed value of the power law index at large scales. Moreover, it is possible to derive 
an analytical form that fits perfectly well the exact numerical derivation. This approximation 
is based on the approximate form (14) for the function Qg. One can then show that Qf (r) = 
t(1 + r 2 /9// 2 (fi, A)) 1 / 2 - r 2 /3//(fi, A). The calculation of the expression (57) then leads to the 
expression, 



2Xai 



d6, (60) 



with 

K = 1+ W)' A = 1 "^A)- (61) 



4.4. Comparison with numerical simulations 

To obtain a smoothed density field from a discrete sample of points is a simple task. To obtain 
the smoothed velocity divergence is not so straightforward and it is worth noting that the smoothing 
is supposed to be done with a volume-weighted scheme, that is the velocities, at any point of the 
considered volume, should be equally weighted, regardless of the local density. 

I used the CDM numerical simulation once again to check the shape of the expected 
distribution function of 6. A grid of 40 3 cubic cells (of size 5 /i _1 Mpc) have been used to define the 
velocity field in the whole simulation, the velocity associated with each cell being the average of the 
velocities of the points contained in the cell. A spherical box of a given radius has been centered on 
each cell. The smoothed velocity in each box is obtained by taking the average of the velocities in the 
cells that are part of the box (weighted by the fraction of the volume of the cell actually contained in 
the box). This procedure insures that the filtering has been made with a volume weighted scheme, 
as soon as the ratio of the radius of the box by the size of the cells is large enough. I made the 
measurement only for a soothing radius of 10 /i _1 Mpc and 15 /i _1 Mpc to keep a high enough ratio 
(and only for a/ao = l). The divergence has been calculated with finite differences using the initial 
grid. 

The resulting probability distribution functions of the velocity divergence are presented in fig. 
6. They appear to be extremely noisy in the tails but the central bump is perfectly well described 
by the analytical form (60) with £1 = 1. The values of erg to be used in (60) can be determined 
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Fig 6. The one— point probability distribution function of the velocity divergence obtained in an Q = 1 CDM 
numerical simulation. The smoothing radius are 15 h Mpc (top panel) and 10 h Mpc (bottom panel) and <7g is 
respectively 0.32 and 0.45. The solid lines correspond to the formula (60) for Q = 1, the dashed lines for Q = 0.3 
and the long dashed lines for Q = 3. 

by suppressing the noisy part of the tails. They are respectively <jg = 0.45 and <jg = 0.32. The 
differences between a and <jg (see table 2) for the same smoothing radius are thought to be due 
to the smoothing scheme used for the velocity field. Indeed, it tends to smooth the velocity field 
in an effective radius greater than the size of the box (especially because of the finite differences), 
and thus to reduce the value of <jg. Moreover the noise is particularly important in the underdense 
regions (for which the velocity is poorly known) and give birth to noisy data in both the low and 
the high density tails (as can be seen with an analysis of the 6—6 correlations). For instance for 
R = 10 /i _1 Mpc about 2% of the measured divergences are in the tails (defined by 6 < —2 or 
6 > 1) while 0.6% of the cells of the grid were completely empty (in such a case the velocity of the 
cell was taken as the average of the velocities in the surrounding cells). Note also that when the 
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tails are suppressed the values found for T3 and T4 are respectively — 1.8 ± 0.15 and 4.9 ± 1.5 in good 
agreement with the theoretical predictions for n = — 1 and £1 = 1 (—1.7 and 4.6). 

The form (60) proves also to be extremely accurate to give the shape of the distribution 
function of the velocity divergence obtained by Juszkiewicz et al. (1993b). 

When the data are compared to the expected shape of the probability distribution function of 
6 for other values of £1, a clear discrepancy can be seen, especially for low values of SI (short dashed 
lines in Fig. 6). Such a function could then be of crucial interest to constrain observationally the 
value of SI. 

6. Conclusion 

In this paper, I have proposed a method to calculate the cumulants of the large-scale 
distribution function of the density and the velocity divergence when a top hat smoothing window 
function is applied to the fields. The generating functions of the cumulants (defined in Eq. [10]) 
are given by the equations (41) and (54) as a function of the shape of the power spectrum. The 
probability distribution functions can then be obtained from the form (18) used with the appropriate 
generating function (p(y). This form has been derived assuming that the ratios {& p ) c / {b 2 ) P ( or 

(@ p ) J ' (0 2 y 1 ) are well approximated by their large-scale limit. Various approximations for the 
probability distribution functions have then been given (Eqs. [47, 50, 52] and [57-61]). The n = —1 
case received a particular focus: it has been shown that the lognormal distribution was particularly 
close to the exact derivation of the density one-point distribution function (although this distribution 
has none of the general properties of the exact results, see text) and a simple analytical fit has also 
been given in that case for the velocity divergence distribution function (Eqs. [60-61]). 

For fields smoothed by a Gaussian window function I proposed an empirical model that 
consists in replacing the index of the power spectrum by an effective one that reproduces correctly 
the skewness. The distribution functions can then be calculated by the results obtained for a top 
hat window function but using the effective index. This is not, however, an exact result. 

Comparisons with numerical simulations proved that these results are extremely accurate in 
the whole quasilinear regime. It implies that the hypothesis of constant ratios, {& p ) c / {& 2 ) P ; is 
extremely powerful in such a regime. Results of numerical simulations (Bouchet & Hernquist 1992, 
Lucchin et al. 1993) indicate that such an hypothesis is no more exact in the fully nonlinear regime 
when a exceeds unity. It may be, however, interesting to know what features are induced by the 
results of the large scale cumulant behaviors when they are applied in the nonlinear regime. This is 
left for a coming paper. 

Applications to observational data are also left for another paper (Kofman & Bernardeau, in 
preparation). 

Note: The Fortran program used to integrate Equation (18) is available upon request (email address: 
fbernard@amoco.saclay.cea.fr). 
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discussion and H. Couchman for providing the numerical simulation. Part of this work was done 
while the author was at the Service de Physique Theorique, CE de Saclay, France 
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